clear; clc; close all;
%% 参数设置
tMax = 200000; % 模拟时间步数
lx=70; ly=100;
nx=100;ny=100;
Player=0; %药物层数
Dspeed=0.5; %扩散速度
x=rand(nx,ny)*lx;
y=rand(nx,ny)*ly/10+1.2*ly;
C=zeros(ly,2);
numParticles = nx*ny; % 颗粒数量
ux=zeros(size(x,1),size(x,2));
uy=zeros(size(x,1),size(x,2));
for tStep=1:tMax
yp=y; %备份y
k=y<ly;
ux=Dspeed*(rand(nx,ny)-0.5).*k+(ux+0.03*rand(nx,ny)).*(1-k);
uy=Dspeed*(rand(nx,ny)-0.5).*k+(uy-0.0001*rand(nx,ny)).*(1-k);
x=x+ux;%+0.5+0.5*sin(tStep/10000*2*pi);
y=y+uy;
x=x.*(x<=lx).*(x>=0)+(x-lx).*(x>lx)+(x+lx).*(x<=0); %循环边界
ku2d=(y<ly).*(yp>ly);
kkup=(y>ly).*(yp<ly); %从下往上越界
kkdown=(y<0).*(yp>0); %从上往下越界
y=(ly-0.2).*ku2d+y.*(1-ku2d);%处理从上往下越界
y=(2*ly-y).*kkup+y.*(1-kkup);%处理从下往上越界
y=(-y).*kkdown+y.*(1-kkdown);%处理从上往下越界
if tStep<2000
Niteval=20;
else
Niteval=200;
end
if mod(tStep,Niteval)==0
clf;clc; tStep
subplot(1,2,1);
hold on
patch([0 lx lx 0],[2*ly/3 2*ly/3,ly ly],[0.4 1 0.3],'EdgeColor', 'none')
patch([0 lx lx 0],[0 0 2*ly/3 2*ly/3],[0.4 0.5 0.8],'EdgeColor', 'none')
% scatter(x,y,'Marker','o','Color','red','MarkerFaceColor','red');
plot(x,y,'r.');
axis equal; axis off
axis([0 lx,0 1.3*ly])
set(gcf,'Color',[1 1 1])
hold off
subplot(1,2,2);
for i=ly:-1:1
C(i,1)=i;
C(i,2)=sum(sum((y>i-1).*(y<i)));
end
barh(C(:,1)/ly,C(:,2));
% axis([0 1000,0,4000])
% axis equal
box on
drawnow
end
end